## find equilibria for economy with 1. no friends, 2) nonassortative matching 3) assortative matching
lowest_quantile<-0.00
highest_quantile<-1.00
low_quantile<-0.25
high_quantile<-0.75
median_quantile<-0.5

integer_quantile<-function(x, quant) 
{
  n<-length(x)
  
  which_one<- min(max(round(quant*n), 1), n)

  x[order(x)][which_one]
}

attach(parameters_new)
mu_vec_tmp<-solution$mu_study_model
gamma_mu_vec_tmp<-gamma_mu(mu_vec_tmp)

mu_lowest_gamma_mu<-mu_vec_tmp[which(gamma_mu_vec_tmp == integer_quantile(gamma_mu_vec_tmp, lowest_quantile))]
mu_low_gamma_mu<-mu_vec_tmp[which(gamma_mu_vec_tmp == integer_quantile(gamma_mu_vec_tmp, low_quantile))]
mu_highest_gamma_mu<-mu_vec_tmp[which(gamma_mu_vec_tmp == integer_quantile(gamma_mu_vec_tmp, highest_quantile))]
mu_high_gamma_mu<-mu_vec_tmp[which(gamma_mu_vec_tmp == integer_quantile(gamma_mu_vec_tmp, high_quantile))]
mu_median_gamma_mu<-mu_vec_tmp[which(gamma_mu_vec_tmp == integer_quantile(gamma_mu_vec_tmp, median_quantile))]

sm_A_match<-matrix(c(0,1,1,0), byrow=T, ncol=2)



s_vec_length<-100
s_not_i_points<-seq(0,7,length=s_vec_length)

lowest_match_equil<-calc_equilibrium(sm_A_match)(rep(mu_lowest_gamma_mu,2))
lowest_friends_BR<-cbind(psi(mu_lowest_gamma_mu)(s_not_i_points), psi(mu_lowest_gamma_mu)(s_not_i_points))

low_match_equil<-calc_equilibrium(sm_A_match)(rep(mu_low_gamma_mu,2))
low_friends_BR<-cbind(psi(mu_low_gamma_mu)(s_not_i_points), psi(mu_low_gamma_mu)(s_not_i_points))

highest_match_equil<-calc_equilibrium(sm_A_match)(rep(mu_highest_gamma_mu,2))
highest_friends_BR<-cbind(psi(mu_highest_gamma_mu)(s_not_i_points), psi(mu_highest_gamma_mu)(s_not_i_points))

high_match_equil<-calc_equilibrium(sm_A_match)(rep(mu_high_gamma_mu,2))
high_friends_BR<-cbind(psi(mu_high_gamma_mu)(s_not_i_points), psi(mu_high_gamma_mu)(s_not_i_points))

median_match_equil<-calc_equilibrium(sm_A_match)(rep(mu_median_gamma_mu,2))
median_friends_BR<-cbind(psi(mu_median_gamma_mu)(s_not_i_points), psi(mu_median_gamma_mu)(s_not_i_points))

plot_range<-range(c(0, as.numeric(lowest_friends_BR), as.numeric(highest_friends_BR)))

## table with policy functions for students with different study types
quantile(solution$mu_study_model, low_quantile)

quantile(solution$mu_study_model,0)

calc_policy_expression<-function(mu_s_in) 
  {
  icept_out <- -theta_vec[3] - theta_vec[4]/gamma_mu(mu_s_in)
  slope_out <- (beta_vec[2]-theta_vec[1] - theta_vec[2]/gamma_mu(mu_s_in))
  
  list(intercept=icept_out, slope=slope_out)
}

## TABLE 3
gamma_mu_s_policy_functions_df<-cbind(
  gamma_mu_s_min=calc_policy_expression(mu_lowest_gamma_mu), 
  gamma_mu_s_25th=calc_policy_expression(mu_low_gamma_mu), 
  gamma_mu_s_50th=calc_policy_expression(mu_median_gamma_mu), 
  gamma_mu_s_75th=calc_policy_expression(mu_high_gamma_mu), 
  gamma_mu_s_max=calc_policy_expression(mu_highest_gamma_mu))
print(xtable(gamma_mu_s_policy_functions_df, caption="tab:policy_functions_by_gamma_mu_s"), file="./output/policy_functions_by_gamma_mu_s.tex")

detach(parameters_new)
